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Bond Orientational Order Parameters in the Crystalline Phases of 
the Classical Yukawa-Wigner Bilayers. 
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PACS 52.27.Lw - Dusty or complex plasmas; plasma crystals 

PACS 82 . 70 . Dd - Colloids 

PACS 64 . 70 . K- - Solid - solid transitions 

Abstract. - We present a study of the structural properties of the crystalline phases for a planar 
bilayer of particles interacting via repulsive Yukawa potentials in the weak screening region. The 
study is done with Monte Carlo computations and the long ranged contributions to energy are 
taken into account with the Ewald method for quasi-two dimensional systems. Two first order 
phase transitions (fluid-solid and solid-solid) and one second order transition (solid-solid) are found 
when the surface density is varied at constant temperature. A particular attention is pay to the 
characteristics of the crystalline phases by the analysis of bond orientational order parameters 
and center-to-center correlations functions. 
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Introduction - Two dimensional crystals of charged 
particles were first observed on a monolayer of electron ad- 
sorbed at the surface of superfluid helium [1] ; multilayers 
df charged ions in laser-beam-cooled trapped plasma [2], 
dusty plasma [3-7] and in colloidal suspensions [8,9] have 
been observed as well. In particular, bilayer systems have 
benefited from a large number of theoretical [11, 13, 14] 
dnd numerical [15-17] studies. 

In colloidal and dusty plasma systems, the interactions 
between particles are quite well approximated by Yukawa 
potentials of the form 



V{r) 



(1) 



where Q is the charge carried by the particles, r the 
distance between particles and \d the Debye screening 
length. In ref. [14], the phase diagram of the Yukawa 
bilayer at T = has been computed and six different 
crystal phases have been identified. Previous numerical 
simulations on bilayers systems have adressed either the 
unscreened Coulomb interaction (A^ — s- oo) [16,17] or the 
strong screening region [16]. 

The crystalline phase diagram of Yukawa bilayer systems 
is quite complicated and several arrangements have been 
found to be thermodynamically stable [2,5,9-17]. Such 
complex phase diagram is interesting from an experimen- 



tal point of view, in colloidal and plasma physics, and from 
theoretical points of view for studies of two-dimensional 
melting [18]. 

In this letter, by using Metropolis Monte Carlo simula- 
tions, we study the structure and properties of the crys- 
talline phases of a bilayer of point particles interacting via 
a repulsive Yukawa potential with large Ad (weak screen- 
ing region). 

Methods and model - The system consists of point 
particles evenly distributed in two parallel layers separated 
by a distance h. Since this work adrcsses the properties 
of the crystalline phase, to avoid irrelevant finite size ef- 
fects due to the shape of the simulation box and to permit 
transitions between crystalline phases, it is necessary to 
allow the shape of the box to change. The Monte Carlo 
simulations are performed at constant N ^ T, h and A, 
where A = L,j:Ly sin^ is the surface of the oblique sim- 
ulation box where the angle 7 is variable ; and Ly 
evolve such as the area A remains constant. A trial move 
of the shape of the simulation box is done every MC-cycle. 
The density of particles in each layer is p = N/2A ; the 
Wigner-Seitz radius a is defined by irpa^ = 1. To fulfill 
electroneutrality and to achieve consistency between the 
Coulomb and Yukawa one component plasma each layer 
carries a neutralizing background with a surface charge 
density a = —NQ/2A (see ref. [19] for more details). In 
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Fig. 1: Bond orientatioiial order parameters for crystalline 
phase IVA as function of a. 

plasma physics, such anisotropic neutraUzing backgrounds 
are responsible of sheath confining potentials that are of- 
ten taken as parabolic potentials [6,20]. 
In all computations presented in this letter h = 1.0, 
Q = 14 and A^i = 10 /i [21] ; for a such large value of the 
Dcbyc screening length Ad, the Yukawa interaction is long 
ranged and the energy of the system must be computed 
with the Ewald method as defined in ref. [19] ; direct trun- 
cation may lead to severe bias of the MC sampling [19,22] . 
The total energy of the bilayer of Yukawa particles in the 
background defined above is given by 

^ . (2) 

where Eb is the sum of the energy of the particle- 
background interactions and the background self energy. 
In Eq.(2), the energy is split into intralayer -Emtro and 
interlayer Einter energies. The number of particles in the 
bilayer is TV = 2048 (1024 particles per layer) and, after 
equilibration, averages are accumulated for 5 x 10^ - 2 x 10^ 
trial moves per particle. A few computations have been 
done with N = 512 and 1058. 

The intralayer gn and interlayer gi2 center to center cor- 
relation functions are defined by 

+ E E (3) 

ieL2 jeL2,j^i 

" leLi JGL2 

Bond orientational order parameters - For a two 

dimensional lattice, with Nb the number of nearest neigh- 



Table 1: Bond orientational order parameters for perfect crys- 
tal phases III, IVA and V. We denote by bi and 62 the 
two primitive vectors of the two dimensional primitive cell. 
For all three phases, | 61 | = | 62 |= bo ; Phase IVA corre- 
sponds to two staggered two dimensional rhombic lattices with 
61.62 = bgcosa, Phase III to staggered square lattices with 
Q = 7r/2 and Phase V to staggered hexagonal lattices with 
a — 7r/3. A'^i, is the number of nearest neighbors in each perfect 
crystal phases and 'I'^ are the values of the bond orientational 
parameters for each perfect crystal phases. 
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bors defined by Voronoi constructions [23], the bond ori- 
entational order parameters are defined by 

*" = I^^E'^^p^*"^^)! (4) 

where 9j is the angle between the bond vector, made by 
a vertex of the lattice and its j*'' nearest neigbor, and an 
arbitrary direction. In Table 1, we give the values of the 
bond orientational order parameter of the two dimensional 
crystal phases (III, IVA and V) for n = 4, 6, 8 and 12 ; by 
symmetry of the crystal phases, VP^ is null for all n odd. 
In the present work, only the three crystalline phases III, 
IVA and V with ] bi ] = ] ^2 | have been found to be stable 
; the rectangular phases I and II for which ] bi [7^] ^2 | is 
found in the phase diagram at T = [14] are replaced by 
a fluid phase at the temperature considered in the present 
work (see below). On Fig.l, we represent as function 
of a. One should note that for the perfect square crystal 
phase III, the Voronoi construction is degenerate : more 
than three Voronoi cells share a given vertex. For phase 
III, the Voronoi construction gives Nb — 4: and the values 
of have to be computed with the four neighbors. 
As shown on Table 1 and Fig.l, the value of 5'° (7r/2; III) 
for n = 4, 6, 12 are not equal to ^°{a ^ 7r/2; IVA) = 1/3, 
because the number of nearest neighbors changes from 
Nb = 4, for a = 7r/2 in phase III, to iV^ = 6 for a ^ 7r/2 
in phase IVA. Therefore, the bond orientational order pa- 
rameter ^"4 is not a convenient choice for the descrip- 
tion of phases III and IVA. However, ^Pg is continuous for 
III^IVA and one has «'g(a -> 7r/2; IVA) = «'g(7r/2; III) = 
1. 



Crystalline Phases of the Classical Yukawa- Wigner Bilayers. 



Table 2: Simulation results for several densities for A'^ = 2048, Q = 14 and Xd = 10 h (h = 1.0). Phases refer to the conventional 
labelling of crystalline phases developed for plasma bilayers (see for instance refs. [2,7, 11]). p is the surface density in each 
layer ; j3U/N is the average energy per particle ; j3 < Einter > /N is the average interlayer energy per particle ; < > and Xn 
axe respectively the average bond orientational order parameters and the susceptibilities given by Eqs.(5,6) and Qg is the angle 
between the two basic primitive vectors of the two-dimensional lattices, computed from the location of peaks in the intralayer 
correlation functions gii(s) (see text and Fig. 5). 
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Fig. 3: Susceptibilities as functions of the surface density for 
N = 2048, Xd = 10 h and Q = 14. For vrp = 1.713, we have 
Xi2 ~ 32.0±4.0. Probability distributions for at vrp — 1.25, 
1.5, 1.713 and 2.0 are given in the inset. 



The average values of bond orientational order parameters 
in the MC samplings are defined by 

1 1 ^"'^ 

<*">=]^(lE]^E6^p(*"%)|) (5) 

where 9ij is the angle between the interparticle vector 
and an arbitrary fixed direction. The number of nearest 
neighbors N^^i of particle i is determined by a Voronoi 
construction. The bond-orientational susceptibilities are 
defined as 

Xn = N{< ^l>-< *„ >2) (6) 

Results and discussion - On Fig. 2, we show the 
bond orientational order parameters and on Fig. 3, the 
susceptibilities for n ~ 4, 6, 8 and 12 ; several results 
for the different thcrmodynamical states are reported in 
TABLE 2. The 'discontinuities' of < ^1/4 > and < "i^s > 
at np ~ 0.2 and of < >, < *8 > and < > at 
np ~ 1.7 are characteristic of discontinuous first order 
phase transitions (the susceptibilities X4: Xe, Xs and 
X12 exhibit also discontinuities at these densities). The 
transition at irp ~ 0.2 corresponds to the transition 
between a disordered fluid phase and the crystalline phase 
III, while the transition at irp ~ 1.7 is the transition 
between crystalline phases IVA and V (see Table 2). 
As outlined before, a perfect ordered square lattice is 
degenerate for the Voronoi construction, therefore, even 
very small fluctuations in position of particles will reduce 
greatly the number of Voronoi cells with four sides {i.e. 
at T 7^ 0). Thus, the values of < ^4 > in phase III, 
computed with Eq.(5), never exceed 1/3, as explained and 
shown on Fig.l ; in the phase III, due to the fluctuations, 
each particle has on average more than four nearest 
neighbors. On Fig. 4(a), we give a snapshot of the bilayer 
for np = 0.75 ; for this configuration, the instantaneous 
bond orientational parameter are : ^'4 ~ 0.32, — 0.01, 
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'I's = 0.74 and ^"12 = 0.11. Since a lattice with a square 
symmetry obviously has also a eight fold symmetry, the 
values of < > are quite large : a moderate value of 
< ^'4 > and a high value of < > are representative 
of phase III (see Tables 1-2 and Figs. 1-2). On Fig.4(b), 
we give a snapshot of the bilayer for np — 1.35 ; for 
this configuration, the instantaneous bond orientational 
parameter are : ^'4 = 0.23, *6 = 0.32, 'ifs = 0.37 and 
*i2 = 0.05. 

The increase of xe and xs at irp ~ 1.35, with no dis- 
continuity for < > and < ^'8 > is associated with a 
continuous second order phase transition. This second 
order phase transition is the transition between the 
crystalline phases III and IVA. 

From the values of the order parameters reported on 
Fig. 2, the bilayer is in a fluid phase for np < 0.2, 
in a staggered square crystalline state (phase III) for 
0.25 < irp < 1.35 ; in a rhombic crystal phase (phase IVA) 
for 1.35 < 7rp < 1.7 and in a triangular-hexagonal state 
(phase V) for irp > 1.7. The densities irp ~ 0.2 - 0.225 
and TTp ~ 1.713 are in the coexistence regions of the 
first order phase transitions. For np > 0.25, this phase 
diagram agrees very well with the one computed in 
ref. [14] for T — ; however, at small density, np < 0.2, 
the crystalline phases I and II, found in ref. [14] are 
replaced by a disordered fluid phase. 
It is interesting to note that < > varies smoothly 
with the density over all the crystalline phase IVA 
(1.35 < wp < 1.7). As shown below, in the analysis of the 
correlation functions, this smooth variation of the order 
parameters is associated with a smooth variation of the 
shape of the primitive cell for phase IVA ; this is also in 
agreement with the analysis performed in ref. [14]. 
In ref. [16], a melting criterion based on the values of the 
bond orientational order parameters < 5*4 > or < > 
(noted Gg in [16]) is proposed : it is argued that melting 
occurs when < ^4 > or < ^Pg >~ 0.45. This criterion 
is clearly incorrect to describe the thermodynamical 
stability of phases III and IVA (cf. Fig.2 and Fig. 4). 
Finite size effects (288 < N < 780), poor sampling 
(2 X 10'^ trial moves per particle) and simulation box with 
fixed shape are certainly responsible for the quite large 
value of Gg observed in ref. [16] in comparison with the 
computation done for perfect crystal phases III and IVA 
reported in Table 1 and Fig.l. 

A few MC histograms P(^'g), normalized to 1, are 
given in inset of Fig. 3 for np — 1.25, 1.5, 1.713 and 
2.0. For np = 1.713, the histogram is double peaked 
since the system is in the coexistence region for this 
density. Single peaked histograms are quite well fitted by 

P(«'n) - Po,n(p)cxp(-iV(«'„- < *„ >)V2Xn). 

On Fig. 5, we show the intralayer 511 (s) and interlayer 
1712(3) center-to-center correlation functions for the three 
crystalline phases III, IVA and V. From the values of 
ffii(s) we can determine the parameters and the structure 
of the primitive cell of the two dimensional Bravais 
lattice in each layer. In crystalline phases, the particles 




Fig. 4: (color online) Snapshots of bilayer systems with Voronoi 
construction for one layer. A'' = 2048. Particles belonging 
to different layer are represented by open and solid circles. 
Voronoi cells with four sides are represented in yellow, those 
with five sides are represented in green, those with six sides 
in white, those with seven sides in red and those with eight 
sides in blue. The sides of the simulation box are represented 
by thick black lines, periodic boundary conditions are applied 
; (a) np = 0.75 (7 = 89.5°) and (b) np = 1.35 (7 = 88.9°). 



fluctuate around their equilibrium position defined by the 
sites of the two dimensional lattice, thus 511 (s) may be 
accurately fitted by a sum of gaussian functions of the 
form gn{s) = go„exp(-(s - S'„)^/2ct^) where Sn is the 
location of the n*'' peak. For the rhombic primitive cell 
(phase IVA), if we denote by 61 and 62 the two primitive 
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Fig. 5: (color online) Intralayer gii(s) and interlayer gi2{s) 
center-to-center correlation functions for the three crystalline 
phases. The correlation functions in black correspond to crys- 
talline phase III (staggered square lattices) ; in blue to phase 
IVA (staggered rhombic lattices) and in red to phase V (stag- 
gered hexagonal lattices). 



vectors of the two dimensional primitive cell and by a the 
angle between these two primitive vectors, then the first 
three peaks are located at Si = bo, S2 ~ 6o-\/2(l — cosa^) 
and S3 = 6o\/2(l + cosa^) with 60 =| bi | = | ^2 |- For 
densities 1.25 < up < 1.6, a fit of the first three peaks of 
511 (s) by a sum of three independent gaussian functions 
gives 60 = 1-80 a and 90° < < 60° (some values of ag 
are reported on Table 2). For 1.35 <TTp< 1.45, < \&i2 > 
is almost null while < ^'g > and < > are almost 
equal, according to the computations reported on Fig.l 
the angle of the rhombic primitive cell is a ~ 80° and 
agree well with the value of ag computed from correlation 
functions (cf. Table 2). 

For phases III and V, a similar analysis of gii(s) can be 
done. For phase III, we found 60 = 1.80 a and ag = 90° 
(the peaks located at S2 and S3 in phase IVA merge into 
a single peak - the second one - in phase III, see Fig. 5(a)) 
; for phase V, wc have bo — 1.90 a and ag = 60° (the 
peaks located at Si and 6*2 in phase IVA merge into a 
single peak). 

The analysis of the peaks of 1712(5) can be done in a similar 
way. For the phase III (staggered square lattices), the 
first two peaks of (712(3) must be located at S'l = bo/^/2 
and S'2 = bo \/5/2. For irp < 1.3, these values are well 
reproduced. For instance, for np = 1.0, fitting the peaks 
by sums of gaussian functions, from 511 (s) we found 
60 = 1-77 a and from 512(3) we have S[ = 1.24 a and 




Fig. 6: (color online) Same as Fig. 4 with Tip = 1.4 (a) A'^ = 
512 and (b) N = 2048. The rhombic shape of the simulation 
box results from the MC trial moves of the shape of the box 
: (a) 7 = 73.8° (< 7 >= 74.4° ± 1.1°) and (b) 7 = 77.4° 
(< 7 >= 76.6° ±0.6°). Whereas < 7 > is strongly correlated 
to the crystalline order of the phase, it may and can not be 
used as an order parameter because of defects (see for instance 
Fig.4 and TABLE 2). 



S'2 = 2.80 a. For np > 1.35, the system is in phase IVA 
and the first two peaks of 1712(5), observed in phase III, 
separate each one into two peaks (see the blue curve in 
Fig.5(b)). 

On Fig. 6, we show snapshots of instantaneous configura- 
tions of the bilayer system with the Voronoi construction 
for one layer. On these figures, Voronoi cells with different 
number of sides are represented with different colors. 
The irregular hexagonal Voronoi cells are typical of the 
rhombic primitive cell of phase IVA. As may be seen on 
these figures, particles belonging to one layer are mainly 
located on the edge of the Voronoi cells of the particles 
belonging to the other layer. 

As it appears on Fig. 6, the main influence of the finite 
size effects is to reduce the number of dislocations per 
surface area. This effect improves the stability of the 
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crystalline phases [18]. It is also worthwhile to note that, 
for the temperature considered in this letter, dislocations 
with five-fold and seven-fold symmetries are grouped 
together in 'neutral' clusters (pairs, quartets, etc.) ; 
larger clusters are observed in larger systems. Similar 
results have been observed experimentally on monolayers 
of strongly coupled dusty plasma monolayers [24] . 
All computations presented here are done in the weak 
screening region of Yukawa potentials where the in- 
teraction is long ranged and the use of Ewald sums 
necessary [19, 22]. The phase diagram shows three 
different phase transitions : two first order transitions, 
at Tvp ~ 0.2 (fluids III) and at Trp ~ 1.7 (IVA^ V) 
and a second order transition at ~ 1.35 (III-^ IVA). 
This phase diagram, except for Trp < 0.2, is in very 
good agreement both qualitatively and quantitatively 
(in the numerical values of the densities of the phase 
coexistence [21]) with the phase diagram computed in 
ref. [14] at T = 0. 
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